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PREFACE 


Li 


This  study  was  initiaLly  motivated  by  Jon  Mike  Smith's  seminar 
"Recent  Developments  in  Modern  Numerical  Metliods  and  Cost  Saving  Te.cli- 
niques"  which  was  held  at  Redstone  Arsenal  on  12,  13,  and  14  February 
1976  [1].  It  was  further  motivated  by  Professor  Charles  A.  Halijak's 
briefing  on  "Numerical  Transforms  and  Digital  Simulation  of  Dynamical 
Systems"  held  during  July  and  August  1976. 

This  report  documents  some  of  the  results  of  the  study  undertaken 
as  a consequence  of  the  previously  mentioned  seminar  and  briefing. 

Though  not  intended  to  be  tutorial,  it  is  hoped  that  it  will  be 
comprehensible  to  those  acquainted  with  z-transforms . 
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SYMBOLS 

Eta,  phase  compensation  or  interpolation  factor 


X( 

) 

Laplace  transform  operator 

n 

Index  or  step  number 

T 

Sample  time  or  time  step 

x("> 

n-th  derivative  of  x(t) 

x(- 

-n) 

n-th  integral  of  x(t) 

2 

delayer,  e 

2( 

) 

z-transform  operator  (Appendix  A) 

I. 


INTRODUCTION 


The  rapid  development  of  digital  computers  since  World  War  II 
has  prompted  a reexamination  of  numerical  computation  techniques  from 
tlie  sampled-data  point  of  view.  A very  interesting  outcome  of  this 
point  of  view  is  that  some  of  the  classical  integrators  are  actually 
phase  shifted  variations  of  the  same  integrator  [1,2].  For  example, 
the  integration  recurrence  relation 


X 


n+1 


X + T [r,x  + (1  - ’,)i  1 

n ' n+1  n 


(1) 


leads  to  the  following  classical  integration  formulas  (Table  1) . 


TABLE  1.  EQUIVALENCE  OF  TUNED  AND 
CLASSICAL  INTEGRATORS 


Phase  Shift  (-,) 

Classical  Name 

0 

Euler  (left  Riemann  sum) 

1/2 

Trapezoidal 

1 

Rectangular  (right  Riemann  sum) 

3/2 

1 

Adams  second  order 

The  fact  that  the  classical  integrators  listed  in  Table  1,  differ 
only  in  their  phase  is  quite  interesting;  even  more  interesting  is  the 
possiblity  of  phase  shifts  leading  to  integrators  between  the  classi- 
cal ones. 

II.  THE  TUNABLE  INTEGRATOR 

The  digital  computer  is  of  necessity  a discrete  time  system 
and,  as  such,  is  amenable  to  analysis  using  z-transforms  [3,  4,  5], 

It  is  well  known  in  z-transforms  that 

Z[f(s)g(s)]  = y(z)f(z),  y(s)  ^ g(s)  . (2) 

It  follows  that 

(3) 

and  taking  g(s)  as  the  input  and  f(s)  as  the  plant,  y(z)  cannot  be 
determined  without  knowing  the  plant,  f(s). 
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If  a continuous  svstem  is  to  be  simulated  on  a digital  computer, 
approximations  must  be  made.  For  example,  if  g(s)  is  not  known  a 
priori  [Figure  1(a)] , Z[ g(s) f (s) ] must  be  approximated  in  order  to 
develop  a recurrence  relation  for  digital  simulation  of  g(s)f(s). 


Figure  1.  Continuous  system. 


If  the  input,  g(s),  is  sampled  and  then  held  by  a zero  order  hold 
[Figure  1(b)],  the  output  after  the  second  sampler  is 

Z f(s)g(z)]  , (4) 

which  may  be  simplified  to 


Given 

f(s)  = 7 , (6) 

2 

the  modified  z-transform  of  1/s  , substituting  Equation  (6)  into 
Equation  (5)»is  (Appendix  B) 


and  Equation  (5)  becomes 


the  tunable  integration  Equation  (1),  in  z-transforrii  notation.  The 
details  of  the  conversion  to  a recurrence  are  shown  in  the  sample 
problem  section. 
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For  a sufficiently  fine  timestep,  T,  it  is  hoped  that 

y(z)  ■ T[  I + (1  - 'i)zjg(z)  (9) 

and  for  a given  input,  g(s) , Equation  (3)  would  allow  a check  on  the 
approximation.  Equation  (9),  for  f(s) , an  integrator. 


An  attempt  might  be  made  to  extend  this  procedure  further  by  inte- 
grating f(s)  by  "inserting"  another  sample  and  hold. 


f(z) 


and,  after  taking  the  modified  z-transform  of  1/s  , 


(10) 


Z[g(s)f(s)]  = T[q  + (1  - q)z]f(z)g(z) 


(11) 


or 

y(z)  « T[q  + (1  - ",)z]g(z) 


(12) 


for  any  f(s)  . 


The  insertion  of  this  sample  and  hold.  Equation  (10),  is  not  physi- 
cally realizable  because  f(z)  is  not  an  input;  this  attempted  extension 
is  not  the  proper  approximation  as  will  be  seen  in  the  next  section. 


111.  TUNABLE  TRAPEZOIDAL  CONVOLUTION 

If  a substitution  of  the  definition  of  the  z-transform. 
Equation  (A-10) , is  made  into  Equation  (2) , it  becomes 


Y,  [ f (nt)*g(nX)]  z''  = 
n=0 


X y(nT)z'^ 

" 00 

X f(nT)z'' 

.n=0 

_n=0 

(13) 


Upon  substitution  of  the  definition  of  convolution  and  taking  the 
Gauchy  product  of  power  series,  the  following  is  obtained; 

“ ,nT  » n 

X z"Jg(T)f(nT  - T)dT  = X z"  X v(kt)f(nT  - kT)  . (14) 

n=0  0 n=0  k=0 


1 


j 

i 

! 


1 


5 


Equating  coefficients  of  powers  of  z,  the  following  is  obtained: 


nT  n 

/ g(T)f(nI  - t)  dr  = X y(kT)f(nT  - kT) 
0 k=0 


(15) 


The  mean  value  theorem  of  the  integral  calculus  guarantees  that 
there  is  some  7 such  that  [6] , 


nT 
/ g(T 


) f(nT  - T)dr  = nT  g(7^^  nT)f(nT  - 7^  nT),  0 < < 1 ; 

(16) 


of  course. 


nT  n-1  (k+l)T 

/ g(T)  f(nT  - T)dT  = ^ / g(x)f(nT  - 'r)dT 

0 k=0  kT 


(17) 


and  there  is  some  7 such  that 
k 


(k+l)T 

I g(-r)  f(nT  - T)d-r  = T g(kT  + 7,  T)f(nT  - kT  - 7 T)  , 

kT 


0 < 7^ 


< 1 


(18) 


Knowing  7 and  7 exist  is  like  knowing  the  solution  to  a differential 
n K 

equation  exists;  knowledge  of  the  existence  of  a solution  is  not  know- 
ledge of  the  solution.  Of  course,  it  is  wise  to  check  for  existence  of 
a solution  before  seeking  a solution. 


For  a sufficiently  small  time  step,  T,  trapezoidal  quadrature 
might  be  used  to  approximate  Equation  (18), 


Tg(kT  - 7 T)f(nT  - kT  - 7 T) 

K K 

+ g(kT)f(nT  - kT)] 


|[g(kT  + T)f(nT  - kT  - T) 

(19) 
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Thtin,  the  right  hand  side  of  Equation  (16)  is  approximately 


T 

j X ls(kT  + T)f(nT  - kT  - T)  + g(kT)f(nT  - kT)]  (20) 

k=0 


and  changing  the  summing  indices 


I X S(kT)  f(nT  - kT)  +I  X g(l<T)f(nT  - kT)  (21) 

k=l  ^ k=0 

and  combining, 

n 

T X g(kT)f(nT  - kT)  [g.f(nT)  + f g(nT)]  . (22) 

k=0 

It  follows  that , Equations  (22)  and  (15), 
n n 

X v(kT)f(nT  - kT)=^  T X g(nT)f(nT  - kT)  - I[g  f(nT)  + f g(nT)]  . 
k=0  k=0  / u u 

(23) 

For  fg  and  g^  both  zero, 

y(kT)  ^ T g(kT)  (24) 

and  from  the  definition  of  the  z-transform  it  follows  that 

y(z)  - T g(z)  (25) 

and 

Z[g(s)f(s)]  ===  T g(z)f(z)  . (26) 

Unfortunately,  requiring  and  g^  both  to  be  zero  is  unusually  restrictive 

because  the  initial  conditions  on  the  input  must  be  zero  and  many  plants 
of  interest  are  excluded  including  the  single  integrator  (Appendix  B) . 
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Substituting  Equation  (23)  into  the  right  hand  side  ot  Equation  (14), 


I II 

T Z z"  Z g(kT)f(nT  - kT)  - 1/2  f„g(nT)  - 1/2  « l'(nT) 
n=0  k=0  I 


(27) 


T Z Z f(nT  kT) 

n=0  k=0 


n-k  T 


fn  Z g(nT)z" 

^ n=0 


+ gg  Z f(nT)z" 

n=0 


(28) 


T 

Z f(nT)z" 

"00 

Z g(nT)z 

T 

0 

^ n=0 

_ n=0 

u 

fo  Z g(nT)/ 

n=0 


+ gn  Z f(nT)z’^ 

n=0 


(29) 


and,  finally,  from  the  definition  of  the  z-trans form,  and  Equation  (29) 


Z[f(s)g(s)]  T f(z)g(z)  - i [fpgCz)  + gQf(z)]  , (30) 


trapezoidal  convolution  [ 7,8,4,]  . It  follows  that 


y(z)  = T 


1 


0 

f(z 


). 


g(z) 


T 

2 


(31) 


If  time  step,  T,  is  to  be  made  as  large  as  possible  for  reasons 
of  speed  in  real  time  simulation,  or  economy  in  Monte  Carlo  studies. 
Equation  (30)  may  limit  how  large  the  time  step,  T,  may  become.  In 
using  tunable  integration  linear  interpolation  is  being  used  for  Equations 
(18), 


Tg(kT  + 7 T)f(nT  - kT  - 7 T)  Tq  g (kT  + T)f(nT  - kT  - T) 

K K K 

+ T(1  - qj^)g(kT)f(nT  - kT) 


(32) 


Of  course,  there  is  no  guarantee  that  a linear  interpolator  will  find 
Equation  (18)  on  the  interval,  0 £ S 1.  It  may  be  necessary  to 

extrapolate.  Extrapolation  may  not  produce  the  desired  values  of  the 
functions;  the  values  at  the  ends  of  the  interval  could  he  equal. 
Hopefully,  this  would  be  the  exception  rather  than  the  rule. 
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■ ■ 


tiki 


Thfin  Equation  (17)  upon  substitution  of  Equation  (32),  would  be 


nT 


n-1 


/ g('t)  f(nT  - T)dT  - T X 

o k=0  L 


’^k+l  T)f(nX  - kT  - T) 


+ (1  - g(kT)f(nT  - kX) 


n n- 1 

= X ][  g(kX)f(nX  - kX)  + X Y (1  - \+i)  g(kT)f(nX  - kX) 


k=l 


k=0 


n-1 


==  r Y g(kX)f(nX  - kX)  + X Y (\  - \+i)g(kT)f(nX  - kX) 


k=0 


k=l 


- X[q^  g^f(nX)  + (1  - yf^g(nX)] 


Assuming 

"'k  ^ ’^k+l 

Equation  (35)  may  be  written 


T Y g(kT)f(nX  kX)  - X[q  g f(nX)  + (1  - q )f  g(nX)] 
k=0 

and  it  follows  that 

Z[f(s)g(s)]  - Xg(z)f(z)  - X[q  g f(z)  + (1  - q )f  g(z)] 

O O CO  O 

Making  a final  assumption. 


- 3 

O 00 


the  following  is  obtained; 


Z [f(s)g(s)]  - Xg(z)f(z)  -X[qg^f(z)  + (1  - q)f^g(z)] 


(33) 


(34) 


(35) 


(36) 


(37) 


(38) 


(39) 


(40) 
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tunable  trapezoidal  convolution.  The  final  assumption,  Kquation  (39), 
does  not  appear  reasonable,  but  the  conditions  under  which  it  may  be 
justified  will  be  discussed  in  tlie  analysis  of  tuning.  From  liquation  (40) 
it  follows  that 

y(z)  - T j^l  - (1  - q)  g(z)  - Tr,g^  . (41) 

In  Section  II,  a digital-to-analog  converter  was  used  to  recon- 
struct an  analog  representation  of  the  input  [Figure  2(a)]. 


(a) 

9 Is)  9 U) 

~9  (s) 

f Is) 

T 

\ -Zlg  (s)f  (s)l 


lb) 


9 Is)  9 lz) 

D/A/D 

~y  lz) 

f lz) 

~y  lz)  f (z) 

T 

Figure  2.  Analog  representation  of  the  input. 

In  this  section  approximations  for  constructing  a digital  representation 
of  an  analog  input  were  developed  [Figure  2(b)].  The  difference  is  one 
of  point  of  view. 


IV.  ANALYSIS  OF  TUNING 


For  a given  plant,  f(s),  and  input,  g(s),  the  relationship 
given  in  Equation  (3)  is  the  exact  digital  representation  of  the 
analog  input.  The  approximation. 


y(z)  - T 


fi-i  -^1 

L 2 f(z)J 


g(z)  - 7 


(42) 


would  be  used  to  develop  recurrence  relations  for  digital  simulation 
(Section  V contains  the  details).  For  a given  plant,  f(s),  and  input, 
g(s),  tlie  ratio  of  Equations  (42)  and  (3)  would  be  a measure  of  tlie 
accuracy  of  the  approximation. 

For  an  integrator. 


f(s) 


I 

s 


(43) 


y(z) 


f z)  g(z)  " 2 ^O 


(44) 
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For  the  input, 


g(t)  = a e"^" 

(45) 

- s + a ’ 

(46) 

and 


g(z)  = 


-at 


1 - z e 

The  exact  representation  would  be 


, . (1  - e )z 

y(z)  = ^ 

1 - z e 


Equation  (42)  yields 


y(z)  - - (1  + z) 


1 - z e 


-at 


2 


(47) 


(48) 


(49) 


2 


(1  + e'^^)z 


1 


z 


-aT 

e 


(50) 


Dividing  Equation  (50)  by  Equation  (48)  to  obtain  the  ratio  of  the 
approximation  to  the  exact  output,  yields 


which  is  just 
2 2 

For  small  angles 


- , aT  _ _2_  . ^ , 


(51) 


(52) 


(53) 


11 


and  the  ratio  would  be 


1 + 


12  ^ 


(54) 


of  course  "a"  need  not  be  real;  consider 
a = iw  , 

then  Equation  (52)  becomes 


(55) 


2 


ctn 


2 


(56) 


the  ratio  for  a sine  wave  input  to  a trapezoidal  integrator.  For 
small  angles 


ctn 


2 


_L  _ 

ur  “ 6 


(57) 


and  the  ratio.  Equation  (56),  is 

1 . . (58) 

1 ^2 

The  ratio.  Equation  (56),  agrees  with  that  obtained  by  the  trans- 
fer function  approach  [9],  which  also  gives  the  ratio  for  a sine  wave 
input  into  some  other  integrators.  Of  particular  interest  is  Simpson  s 
rule. 


2 

T 1 + 4z  -f  z 
3 (1  + z)(l  - z) 

for  which  the  ratio  is 

uT  / 2 + cos  \ 

3 \ sin  iJT  / 

and  for  small  angles.  Equation  (60)  becomes 


1 


(59) 


(60) 


(61) 


j 

I 

i 


12 


'I 

I 


I 


J 


It  is  noted  that  for  too  large  a time  step,  T,  the  ratio  for  an 
exponential  input,  Equation  (52)  would  be  too  large  and  for  a sine 
wave  input.  Equation  (56)  would  be  too  small.  This  can  be  corrected 
by  tunable  integration.  Using  tunable  convolution.  Equation  (40), 

••  ••  ^ 

^ T[’V(1  - g - J—  S. (62) 

- _aT 


and  the  ratio  of  Equations  (62)  and  (48)  is 


aT 


-aT 

e 

-aT  , 
e - 1 


(63) 


For  "a"  real.  Equation  (63)  becomes 


aT 


cosh  + sinh  — 
aT 

2 sinh  — 


- U 


(64) 


which  simplifies  to 


f [ctnh  f +2(i-  ,)]  . 

Solving  Equation  (65)  for  q , when  the  ratio  is  one, 

11/  aT  2 \ 

q = - + - (ctn  — “ ) . 

For  aT  < 1,  Equation  (64)  becomes  approximately 

,2 


(65) 


(66) 


1 + aT  (^"h)  + ■^^2“  ■ 


and  Equation  (66)  becomes 
1 ^ aT 

^ - 2 ^ l2  • 

For  aT  » 1 , 


(67) 


(68) 


-if 


(69) 


I 


i 

I 


4 
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For  "a"  imaginary,  Equation  (63)  becomes 


/off  \ 

r cJT 

(1  \1 

1 

( 2 ) 

1 ctn  — + 2i. 

(2  ■ 

• 

(70)  1 

Because  Equation  (70)  is  complex,  multiply  by  the  complex  conjugate 
and  take  the  square  root,  which  yields 


(-)|^ctn  - + 4 (2  -Tl)  J 


2-'  I 


(71) 


for  the  amplitude  ratio,  and  dividing  the  imaginary  part  by  the  real, 
yields 


tan"^  [ 2 (I  - q ) tan  ^ ] 


(72) 


the  phase  error  for  a sine  wave  input  to  a tunable  integrator. 
For  zero  phase  error. 


h 2 ’ (73) 

the  trapezoidal-integrator,  and  Equation  (71)  reduces  to  Equation  (56). 


For  zero  amplitude  error,  the  ratio  should  be  one;  therefore, 
setting  Equation  (71)  equal  to  one 


= 1^  ctn  - + 4 ( - - q ) J 


(74) 


and  solving  for  q. 


_l^ir/2\^  21^] 

2 2 [(  uff  ) " 2 J 


(75) 


For  small  angles 


2 - 2 off  / I (JT  6 • • • / J 

. Fi  . 2 

2 ■ L6  144  J 


(76) 


(77) 
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and  for  wT  < 1. 


^h(i) 


that  is , 


0.9082482905 

0.0917517095 


To  determine  the  ratio  and  phase  error  which  occurs  from  ignoring 
the  dependence  of  t)  on  wT,  Equation  (78)  is  substituted  into  Equations 
(71)  and  (72),  and  the  following  are  obtained 


(f)[ 


2 ^ .2 
- + ctn 


[^(i) 


] . 


for  the  amplitude  ratio  and  phase  error,  respectively,  the  ratio  for 
a small  angle  approximation  would  be 

1 * ^ (82) 

which  compares  very  closely  with  Simpson's  rule,  Equation  (61),  but  the 
phase  error  would  be 


Because  w = 


where  "p"  is  the  period,  it  follows  that 

P _ 2n 
T (JI 


is  the  number  of  samples  per  period.  At  the  Shannon  sampling  limit, 


therefore, 
off  = Jt 

The  amplitude  ratio.  Equation  (71),  is  then 


at  the  Shannon  limit.  For  the  trapezoidal  integrator  the  ratio  would 
be  zero;  the  cost  of  zero  phase  error. 

Simpson's  rule  also  has  no  phase  error  but  at  the  Shannon  limit 
the  ratio.  Equation  (60),  would  be  infinity. 

For  Equation  (78),  Equation  (80)  gives  the  ratio  1.28  ...,  a 
28%  amplitude  error  for  small  angle  tuning  at  the  Shannon  limit 
better  than  Simpon's  rule. 

For  unity  amplitude  ratio  at  the  Shannon  limit. 


that  is, 

I 0.8183098861 
I 0.1816901139 

For  Equation  (83),  Equations  (71)  and  (72)  become  respectively. 


(84) 


(85) 


(86) 
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for  small  angles,  the  ratio.  Equation  (86),  would  be  approximately 


1 - 

15.3 

not  much  better  than  trapezoidal,  Equation  (58). 

Thus  far,  only  exponential  inputs  to  an  integrator  have  been 
considered.  Since  functions  may  be  expressed  in  Taylor  series. 


f(t)  = J , 


another  input  of  interest  would  be  powers  of  t. 


S(t)  = 1 , 


g(s)  =7 


and  trapezoidal  convolution  yields 


(i-i) 


(1  - z) 


which  is  exact  and  the  ratio  is  one. 


g(t)  = t , 


g(s)  = 


and  trapezoidal  convolution  yields 

J±  . i\  = 

\ s^  ® / 2(1  - z)^ 


which  is  also  exact  and  again  the  ratio  is  one.  This  in  not  surprising 
because  a trapezoidal  integrator  can  integrate  a constant  or  a ramp 
perfectly  for  any  size  time  step. 


For 


g(t)  = t" 
2! 

g(s)  = ^ 


and 


»)■ 


zd  + 4z  + z^) 
3(1  - z)^ 


and,  trapezoidal  convolution  yields 


B-s)- 


z(l  + 2z  + z^) 

2(1  - 


Dividing  Equation  (98)  by  Equation  (97)  yields 
2 


3 

2 


1 + 2z  + z 
1 + 4z  + z^ 


(95) 

(96) 


(97) 


(98) 


(99) 


The  "z"s  do  not  cancel  out  as  in  the  previous  examples.  Applying 
the  final  value  theorem,  in  the  limit  as  z approaches  one.  Equation  (99) 
approaches  one;  but  applying  the  initial  value  theorem,  as  z approaches 
zero.  Equation  (99)  approaches  three  halves,  a 507„  error. 

Substituting  the  definition  of  the  z-transform  into 
the  following  is  obtained: 


(1  - z) 


00 


z 

n=0 


t = nT 


==  Ttq  + (1  - T))z]  X g(nT)z" 
n=0 


(101) 
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which  mav  bo  reduced  to 


00  —1  / \ / 00  -1/  \ / 

I ^ {^}/ 

t=nT  t=(n-l)T 


n ^ 

ilT  X 8(nT)z  - (1  - ’i)T  X g(nT  - T)z''  . 


(102) 


8(t)  . 


(103) 


y(s)  =-^ 


(104) 


(^)  -T 


(105) 


Equation  (102)  becomes,  after  equating  coefficients  of  like  powers  of  z. 


I [(nT)^  - (nX  - T)^]  ==  | [q(nT)^  + (1  - q)  (nX  - X)^]  , n > 0 . 

(106) 

taking  the  difference  between  the  right  side  and  the  left  side  to  deter- 
mine the  error  per  time  step,  the  following  is  obtained; 


wliich  implies,  for  zero  error,  that 


(107) 


1 / n - 3 

T “ 2 I T 

\n 


(108) 
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For  n = 1,  T)  = — and  it  is  clear  wliy  the  initial  value  for  the  trapezoidal 

integrator,  is  three  halves.  It  is  noted  tliat  , is  different  for 

each  step. 

Summing  the  error  per  step  for  m steps 


gives  the  total  error  for  m steps, 

[(f  ■ ^)f  (™  + 1)  + (f  -i)"]  . 

For  the  total  error  for  m steps  to  be  zero. 


(109) 


m > 0 


(110) 


T = 


1 /m  - 4 


(111) 


n = - (m  + 1)  , 


Equation  (111)  is  in  agreement  with  Equation  (108).  Therefore,  the  n 
which  tunes  for  m steps  is  exact  for  the  middle  step,  and  the  errors  for 
each  step  are  averaged  out.  The  mean  value  theorem  applies  for  one  step. 
Equation  (108),  or  m steps  taken  together.  Equation  (111).  Because  it 
is  not  practical  to  adjust  the  timing  for  each  step,  an  overall  tuning 
using  Equation  (111)  is  indicated.  This  leaves  the  initial  value  of  the 
ratio  uncorrected. 


In  the  earlier  examples  where  the  z's  cancelled  out  when  the  ratio 
was  taken,  it  would  appear  that  n is  the  same  for  all  steps,  because 
there  was  no  dependence  on  the  step  number.  That  this  occurs  for  some 
functions  seems  to  be  the  power  behind  tunable  integration. 

The  analysis  of  an  exponential  input  to  an  Integrator  would  also 
apply  to  a constant  into  a single  pole  filter  because  the  product  of 
the  Laplace  transforms  is  the  same.  A damped  sine  wave  and  a single  pole 
filter,  etc.  could  also  be  analyzed. 


Ik 


V.  SOME  SAMPLE  PROBLEMS 


The  following  relationships  from  Laplace  transforms  will  be 
found  helpful  in  developing  recurrence  relations; 


- "I 


n-1 


I U'“'(t) 


and 


t=0 

n „(-f) 


n-i-1  (i) 

s X 

o 


(112) 


U.  c?  e 


" iTi  ^ 


i+l 


(113) 


When  the  initial  conditions  are  nonzero,  one  should  proceed  from  the 
differential  equation  using  these  relationships,  and  not  directly  from 
the  transfer  function. 


The  z-transform  of 


.nT 


Z g(s)1  = Z f (ll'^) 

° J n=0 


because  a constaiat  in  the  frequency  domain  transforms  to  a Dirac  delta 
in  the  time  domain.  From  the  properties  of  the  Dirac  delta  (Appendix  C) 
Equation  (114)  becomes 


Z [ S(s)]  = X 


(115) 


n=0 


= g(z) 


(116) 


This  relationship,  Eq..ation  (116),  will  be  used  in  incorporating 
initial  conditions  into  the  recurrences. 


A.  Single  Integration 


x(s)  = 


From  Equation  (113)  for  n = 1, 
x(s)  + X 


(117) 


Taking  the  z-transfnrm  of  Equation  (117),  tlie  following  is  obtained; 

x(x,  -z(a^l)  +x_^z(i)  . (1181 


i 
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Using  trapezoidal  convolution,  Equation  (30), 

From  the  definition  of  the  z-transform 
00  00 

(1  - z)  V x(nT)z'^  ^ (1  + z)  ^ x(nT)z‘^  " ^ ^r, 


(119) 


(120) 


^ [x(nT)z'^  - x(nT)z'^^^]  — ^ [x(nT)z'^  + x(nT)z'^^^] 


n=0 


n=0 


T • , 

- X + X 

2 o o 


(121) 


00  00  00  00 

^ x(nT)z”  - ^ x(nT  - T)z"  ^ x(nT)z"  + ? ^ x(nT  - T)z" 
n=0  n=l  n=0  n=l 


(122) 


1 • ^ 

- - X + X 

2 o o 


Equating  coefficients  of  powers  of  z, 

X = X n = 0 

o o 

^n  = Vl  + 2 (^n  ^ Vl^>  " > 

the  trapezoidal  integrator. 

In  a similar  fashion  for  tunable  convolution.  Equation  (40), 

and  Equation  (125)  becomes 

(1  - z)x(z)  - T[  q + (1  - q)z]jl(z)  - Tqx^  + x^ 


(123) 

(124) 


(125) 


(126) 
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It  follows  tliat,  after  manipulations  similar  to  those  in  Equations 
(120)  to  (122), 


X = X 
o o 


X = X + T I >1  X + ( 1 - h )x  , ] , n > 0 , 

n n-1  n n-1 


the  tunable  integrator. 

To  illustrate  the  difference  between  a constant  in  the  time  and 
frequency  domain,  the  following  differential  equation  is  considered 

i(t)  = k . (12^ 

In  the  frequency  domain.  Equation  (129)  becomes 

s x(s)  - X = - (13C 


and  solving  for  x(s). 


k 

- + X 

s o 


Transforming  Equation  (131)  back  to  the  time  domain,  the  following  is 
obtained: 


t 

x(t)  = ^k  + x^  5(t)Ji 


k / dT  + X 


= kt  + X 


Instead,  tlie  z-transform  is  taken  of 


= — + — , 
s 
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which  would  be 


x(z)  = 


(1  - z) 


2 1- 


(136) 


kTz+(l-z)x 

(1  - z)^ 


(137) 


The  recurrence  would  be 


= X 

o o 


(138) 


X,  = X + kT 
1 o 


X = 2x  ,-x  -,n>l 

n n-1  n-2 


(140) 


This  recurrence  may  be  unexpected  but  it  should  be  tried  because  it  works. 
Consider 


1 r kTz  . 1 


(141) 


and  noting  that 


(142) 


one  may  also  have  the  recurrence 


X = X 

o o 


, n = 0 


(143) 


X =x  ,+kT,n>0 
n n-1 


(144) 


This  recurrence  lias  an  advantage  in  that  no  difference  is  taken  and 
only  one  step  is  required  to  start,  though  both  are  proper,  that  is, 
correct  recurrences. 


w 


B.  Double  Integration 

From  Equation  (113)  for  n = 2, 


t V \ X X 

»<.)  . ii|i  . -I  . ^ 

s s 


(155) 


and  taking  the  z-transform 


(156) 


Using  trapezoidal  convolution.  Equation  (30), 

Tz 


\ T Tzx 

2 ) ■ I ^ 

/ (1  - 2) 


and  Equation  (156)  becomes 


(157) 


(1  - z)^  x(z)  — T^z'x"(z)  - ^ z’x’  + Tzx'  + (1  - z)x 

zoo  c 

Equating  coefficients  of  powers  of  z, 

X = X , n = 0 , 

o o ’ 


(158) 


(159) 


X,  = X + T X H — — X , n = 1 , 

1 o o 2 o ’ 


(160) 


and 


2 . 

X =2x  ,-x  „+T  X ,,n>l 

n n-1  n-2  n-1 


(161) 


In  general,  the  number  of  steps  before  the  recurrence  for  x^  may 

be  applied  is  equal  to  the  order  of  the  denominator  of  the  "transfer 
function."  The  required  information  for  the  start-up  steps  is  contained 
in  the  initial  conditions,  even  when  they  are  zero  (Appendix  D). 

The  recurrence,  Equation  (161),  does  not  appear  desirable  numerically 
because  a difference  is  required  and  it  takes  two  steps  to  start. 

Because  no  feedback  is  involved,  two  single  integrators  could  be  used 
to  overcome  these  difficulties; 
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, n = 0 


(162) 


X = X 

o o 


X = X + ^ (x'  + x'  , ) , n > 0 

n n-1  2 n n-1 


X = X 

o o 


, n = 0 


= Vl  + 2 Vl^’  ° • 


C.  A Single  Pole  Filter 


(163) 


(164) 


(165) 


Consider  the  following  differential  equation 


x(t)  = -a  x(t)  + g(t) 


This  would  transform  to 


s x(s)  - x^  = -a  x(s)  + g(s) 


(166) 


(167) 


and  finally 


x(s)  = 


g(s)  + x^ 


(168) 


in  tlie  frequency  domain.  Taking  the  z-transform  of  Equation  (168), 


x(z)  = Z ) + X Z ( — - — ) 

\s+a/  o \s+a/ 

Using  trapezoidal  convolution.  Equation  (30), 


(169) 


.(ffti-Kifit;)-' 


")  ■ 2 


1 - ze 


(170) 


and  Equation  (169)  becomes 


.1  -aT,  , X T , -al\  . x T 

(1  - ze  )x(z)  - - (1  + ze  )g(z)  " 2 


Equating  coefficients  of  powers  of  z, 


X = X 

o o 


^ -aT 
X — e X 
n 


n-l  + f [Sn 


-aT 
+ e g 


, n = 0 

n-l]>  > 0 • 


(171) 


(172) 


(173) 


The  z-transform  of  the  following  equation  might  have  been  used  instead 
g(s)  - a x(s)  + x^^ 


x(s)  = 

which  is 

■•w  = ^ ) - -o  ^ (t)  • 

The  recurrence  would  then  be 

X = X , n = 0 

o o 

X =sx  ,+^Tg  - ax  + g -ax  ],n>0 

n n-l  2 L n n °n-l  n-lj 


and  solving  for  x^. 


X 

n 


i-f 


^*T 


T / fnj_f^l 
aT  / ^n-1  2 \ , . aT 


, n > 0 


It  is  noted  that  flOl 


= (1/1)  Fade  approximation  for  e 


-aT 


(174) 

(175) 

(176) 

(177) 


(178) 


(179) 


' ( 


1 
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and 


j = (O/l)  Pad!  approximation  for  e 


2 


(180) 


Instead  of  solving  for  x , a past  value,  x , might  be  used, 
that  is, 


<—  X i+?  g - ax  ,+g  ,-ax  , 
n n-1  2 L n n-1  '’n-1  n-lj 


which  would  lead  to 


X ( 1 - aT)  X i+^Tg  + g 
n n-1  2 L n ®n-l J 


and  it  is  noted  that  [ 10  ] 


and 


(1  - aT)  = (1/0)  Fade  approximation  for  e 


-aT 


-aT 


1 - (0/0)  Fade  approximation  for  e 
The  declining  accuracy  of  the  approximation  is  also  noted. 


(181) 


(182) 


(183) 


(184) 


Additionally,  Equation  (182)  could  be  used  to  predict  x , substi- 
tute Equation  (182)  into  the  right  side  of  Equation  (177),  " 

% “ %-i  + I [e„  + Vi]  - f [“  - ">Vi  I (8„  + Vi> 


+ ’‘n-l] 
and  simplyfing 


(185) 


( 


1 - aT  + 


(186) 


where 


2 

^1  - aT  + I = (2/0)  Fade  for  e 


and 


(-f)  = 


(1/0)  Fade  for  e 


2 


(187) 


(188) 
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In  these  cases  where  "a"  is  time  varviiig,  and  it  is  not  desired  to 
"ST 

compute  e at  each  time  step,  the  desired  Pade'  approximation  should 
be  substituted  directly  into  Equation  (173). 

D.  A Forced  Damped  Oscillator 

The  differential  equation  is 

x(t)  + 2 C w x(t)  + 10  ^ x(t)  = g(t)  (189) 

o o 

whereto^  is  the  undamped  natural  frequency  and  *•  is  the  damping  factor. 
From  Equation  (112) 

x(s)  = s x(s)  - X (190) 

o 


and 


’x(s)=s^x(s)-sx-x  . (191) 

o o 

Substituting  Equations  (190)  and  (191)  into  Equation  (189)  and  solving 
for  x(s) , 

g(s)  + (s+2!'to)x  + x 

’0^0  O 

x(s)  = 2 2 ' (1-92) 

s -t-2!'(o  s + u 
" o o 

It  is  noted  that 


x(s)  ^ 

g(s) 


s^  + 2 


_1 

( w s + w 
o o 


(193) 


only  if  X =0  and  x =0. 
o o 

Three  possible  ways  of  implementing  Equations  (192)  and  (193)  arc 
sliown  in  Figure  3.  Figure  3(c)  would  be  a proper  implementation  on 
an  analog  computer  and  is  usually  tlie  approach  taken  in  digital  simu- 
lation. 


30 
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A perusal  of  Appendix  B will  reveal  no  z-transform  for  hquation 
(193)  but 

s^  + 2£,u)  s + w^=(s+  - + (194) 

^ o o o 


where 


1/2 

w = (1  - f;  ) '-g 

is  the  damped  natural  frequency.  For  convenience,  let 

a = C U) 

^ o 


(195) 


(196) 


and  for 


2 2 
(s  + a)  + u 

2 2 

0 

S + Ul 

(0,  1) 

2 2 
(s  + a)  + w 

1 

(s  + a)^ 

,2  2 

1 < 

(s  + a)  - w 

Again  referring  to  Appendix  B,  it  will  be  noted  there  are  four 
possible  z-transforms  for  Equation  (193)  and  the  damping  factor  deter- 
mines which  should  be  used. 


The  case  where  there  is  no  damping,  - 0,  and  no  forcing  function, 
g(t)  = 0,  may  be  of  interest  if  a free  running  oscillator  is  required. 

In  this  case.  Equation  (192)  reduces  to 


x(s) 


SX  + X 

O o 


2 

s 


+ w 


(198) 


and  taking  the  z-transform 


x(z)  = x^  z(  2"'— 2)  2 ^ 2) 

\s+tl/  \S+U)/ 


(198) 
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whicli  becomes 


x(z)  = 


(1  - z cos  cJT)  + x^  (z  sin  u/T)/cj 
1 - (2  cos  ciT)z  + z^ 


(200) 


Substituting  the  definition  of  the  z-transform  for  x(z),  Equation  (200) 
becomes 


00 

+ I M 


L x(nT)z  - 2 cos  (JT  ^ x(nT  - T)z  + ^ x(nT  - 2T)z 
n=0  n=l  n=2 

= X + [-(cos  wT)x  + X (sin  aJr)/oj]z 
o o o 


(201) 


which  may  be  written 


n=0  L 


x(nT)  - (2  cos  tJr)x(nT  - T)  + x(nT  - 2T) 


+ (2  cos  tJr)x(-T)  - x(-2T)  - x(-T)z 


= x^  + [-(cos  + x^(sin  oJr)/u]z 


(202) 


Equating  coefficients  of  power  of  z, 


X = X 
o o 


, n = 0 


(203) 


Xj^  = x^  (cos  cJT)  + x^  (sin  tJl')/u),  n = 1 


(204) 


X = (2  cos  (dT)x  , - X - 
n ' n-1  n-2 


, n > 1 


(205) 


The  state  transition  method  yields  the  following  recurrences  [5]  , 


X = (cos  (j.fr)x  , - (tj  sin  cdT)x  , , 

n n-1  n-1  ’ 


(206) 


t _ / .tdT  \ ^ ^ ca')x 

n \ 0)  / n-1  n-1 
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In  comparing  Equation  (205)  with  Equation;  (206)  and  (207),  it  is  noted 
that  the  state  transition  metliod  requires  four  multiplications,  one 
addition,  and  one  subtraction  and  two  stores,  while  Equation  (205) 
requires  one  multiplication,  one  subtraction,  and  two  stores  after  a 
one  step  start  up.  Equation  (204).  Tlie  start  up  step.  Equation  (204), 
is  the  same  as  Equation  (207). 

These  recurrences  are  exact  because  there  was  no  requirement  to 
invoke  any  approximations.  The  exactness  is  easily  demonstrated.  If 

X = cos(iJT  + 0)  (208) 

it  follows  that 

X = -w  sin(wt  + <!>)  (209) 

and  for  t = nT , 

x^  = cos  , n = 0 , (210) 

Xj^  = cos  ufT  cos  <t’  - sin  u/T  sin  , 


= cos  (ufT  + ft) 


, n = 1 


(211) 


and 


x^^j^  = 2 cos  uil  cos(ncJr  + 0)  - cos((n  - l)cJT  + 0)  = cos  ((n  + l)uT 


+ <))),  n > 1 


(212) 


Of  the  four  possibilities,  Equation  (197),  the  case  where  the 
damping  factor,  t,  is  on  the  interval  between  zero  and  one  is  t!ie  most 
useful . 


x(z)  ^ T Z 


2 1 
, (s  + a)  + CO  / 


>8(")  - T 


2 ®o 


+ X Z 
o 


s + 2a 


2 2 
(s  + a)  + CO  / 


+ X Z 
o 


2 2 
, (s  + a)  -t  u , 


(213) 
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taking  the  z- trails  forms,  the  following  is  obtained; 


The  same  procedure  would  be  followed  in  developing  recurrences 
for  the  other  values  of  the  damping  factor,  , in  Equation  (1S17). 

One  problem,  the  values  of  x and  x"  may  be  desired.  This  leads 

n n 

to  a combination  of  Figures  3a  and  3b  (Figure  4). 


Figure  4.  Figures  3(a)  and  3(b)  combined, 
Writing  Equation  (189)  as  follows 

2 

k(t)  + 2 t u)^  x(t)  = g(t)  - x(t) 

and  using  Equation  (112),  the  following  is  obtained; 

2 

(s  + 2 /•  w )x(s)  = g(s)  - u x(s)  + X 


(218) 


and  finally 

g(s)  - w ^ x(s)  + X 

— TT  T3 ^ • ^220) 

" o 

This  is  the  single  pole  filter  discussed  in  the  prev'ious  section. 
Taking  the  z-transform  of  Equation  (220),  the  following  is  obtained; 

(1  - z e‘^'*'^)x(z)  (1  + z e'^'^'^)[g(z)  - x(z)] 

(22n 

-?(g  - w^x)+x 

2 \ o o o / o 

and  the  recurrence  would  be, 


X = X 

o o 


, n = 0 
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and 


. _ -2aT  T r 

n n-1  2 I n 


,2  ^ -2aT 

w X + e 
o n 


(«n-l  - Vl)]>  " 


> 0. 


(223) 


would  be  computed  using  Equations  (216)  and  (217);  then  x^  would  be 


computed  using  Equation  (223)  and  finally  x'^  using  Equation  (224), 


x‘=g  - 2('Wx  - u X 

n n ^ o n on 


(224) 


Because  the  proper  recurrence  for  x^  is  a function  of  the  damping 


factor,  (* , the  preceding  approach  may  be  too  complex.  Let 


T 

X = x +-(x  + x ) 

n n-1  2 n n-1 


(225) 


and  substituting  Equation  (225)  into  Equation  (223),  the  following  is 
obtained; 


-2aT  . , T 

X = e X , + — 

n n-1  2 


e - w 

®n  o 


fx  + ^ (x  + X )1 
L n-1  2 n n-1  J 


(226) 


-2aT  . 2 . 

+ e (g  1 " ^ 1 ) 

'“n-1  o n-1 


and  solving  for  x , 
n 


r 2" 

/u)  T\ 

^n-l  ^ f 

r -2aT  ,2  ^ -2aT,  1 

K + " Vl  - ‘"o  " >"n-lJ 

X = 

n 

'vf 
^ 2 ' 

(227) 

This  would  correspond 

to  Figure 

3(b). 

Another  way  to  develop  a recurrence  would  be  to  integrate  twice; 
this  is  the  more  traditional  approach.  Let 


T . . 

X = X , + — (x‘  + x‘  , ) 
n n-1  2 n n-1 


(228) 
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and 


X = x + — (x  + x ,) 
n n 2 n n-r 


where 


x'=g  - 2/'Wx  - CO  X 

n n > o n on 


(230) 


Because  there  are  three  equations  and  three  unknowns,  x'  , x , and  x , 

^ n n n 

Equation  (230)  is  substituted;  then  Equation  (229)  is  substituted  into 

Equation  (228)  and  solved  for  x^^, 


1 - (^)  . 
^ ^ 


, T 8n  ^ «n-l  ' ^ 

2 ^.-1  ^ 2 n~ 


1 + ‘'OJ  T -f 
o 


m 


(231) 


This  would  correspond  to  Figure  3(c). 

E.  Newtonian  Drag 

A nonlinear  problem  of  interest  in  missile  and  aircraft 
simulations  is  the  Riccati  equation  representing  Newtonian  drag  with 
acceleration.  Equation  (232); 

d(t)  = -k(t)  u^(t)  + g(t)  , (232: 

which  in  the  frequency  domain  may  be  written 


s u(s)  - u^  = - X[k(t)u  (t)]  + g(s) 


- £[ k(t)u  (t)]  + g(s) 


(233) 


i 2 V.) 


Taking  the  z-transform  of  Equation  (23A)  the  following  is  obtaineil; 


■X 

LV 


I 


With  the  help  of  trapezoidal  convolution.  Equation  (30),  Equation  (235) 
becomes 


(1  - z)  u(z)  = - (1  + z)  Z 


+ u 


g(s)  - £[k<t)u^t>]  -f  u/] 

(236) 


Equating  coefficients  of  powers  of  z,  the  following  is  obtained: 


u = u , n = 0 

o o 


(237) 


and 


u = u + g u^  + k u^.l  ,n>0 

n n-1  2 L n n-lj  2 L n n n-1  n-lj  ’ 

(238) 


T,  \ 

— ku  +u  - I1-—  k u iju  ,+r  g + g , 

2 n n n \ 2 n-1  n-1  f n-1  2 n n-1 


Letting 


= n = (^  - 2 ^-1  Vl)  Vl  ■'f  + 8n-l 


(239) 


(240) 


and  solving  for  u^,  the  following  is  obtained: 


2 c 


1+  (1  + 2Tk^  c^) 


n , n > 0 

T72 


(241) 


This  recurrence.  Equation  (241),  seems  a bit  complex.  Letting 


c 

II 

1 

n-1 

+ 

\ n 

- u 

n 

-li 

f 

= -k 

2 

+ 

2k 

u , 

u 

+ k f u - u 

n 

n-1 

n 

n-1 

n 

n\  n n 

; th  at 

(\ 

’ 1 

f "" 

negligible 

-k 

n 

2 

n-1 

+ 

2k 

n 

^ 1 
n-l 

u 

n 

. 

(242) 

(243) 


(244) 
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Substituting  Equation  (244)  into  Equation  (238)  and  solving  for 
yields 


2 (‘^n  ~ '^n-l)  %-l]  2 [^n  ^n-1  ] 


1 + Tk  u , 
n n-1 


, n > 0 (245) 
(245) 


A Small  angle  approximation  to  the  analytic  solution  of  Equation 
(232),  for  g and  k constant,  is  [ llj , 


^o  + 

1 + t k u 


(246) 


It  would  appear  that  Equation  (245)  is  a better  approximation  than 
Equation  (241)  when  compared  with  the  analytic  solution  111]. 


Another  possibility  is  to  use  the  relationship  between  the  Riccati 
equation  and  the  linear  differential  equation  of  second  order  [12]  in 
Equation  (232)  which  becomes  upon  transformation. 


u(t)  = 


k(t)y(t) 


which  becomes 


V‘(t)  = y(t)  + k(t)  g(t)  y(t) 


(247) 


(248) 


The  procedure  would  be  similar  to  that  in  the  treatment  of  tlie  forced 
damped  oscillator,  and  will  not  be  pursued  here. 

Equations  (241)  or  (245)  would  be  of  use  only  in  a one-degree-of- 
freedom  simulation,  say  rising  or  falling.  In  general. 


u = -kVu  + g 


(249) 


V = -kVv  + h 


(250) 


■kVw  + i 


(251) 


where 


V = (u  + V + w ) 


(252) 


40 


In  body-fixnd  coordinates,  u — V, 


1 


“ - - h izf  * (?fj  ' ^ . 

(253) 

j. 

^ {u}  ■"(?)]  ^ ' 

(254) 

1 

f {z)  ^ " • 

(255) 

Let 

J, 

k'  = k(l  + 

(256) 

where 

1 

-^[(i)^(z)]' 

(257) 

t 

is  the  total  angle-of-attack.  Then  Equations  (249), 

(250),  and  (251) 

become 

u = -k  u + g 

(258) 

V = “k'uv  + h 

(259) 

w = -k  *uw  + i 

(260) 

The  recurrence  for  Equation  (260)  would  be  the  same  as  for  Equation  (259) 
that  is,  the  obvious  substitutions  are  made  in  Equation  (261).  The 
rate  equations  would  be  handled  in  a similar  manner. 

VI.  CONCLUSIONS  AND  RECOMMENDATIONS 

Those  who  didn't  know  have  probably  guessed  tiiat  a potpourri 
is  a "disli  of  severall  meates  boyled  and  stued  togctiier."*  Some  might 
be  of  the  opinion  it  should  have  been  allowed  to  cook  a little  longer, 
but  this  report  was  not  intended  as  a final  summary;  it  was  intended 
to  document  tlie  results  of  an  initial  investigation  and  suggest  avenues 
for  further  development.  They  are  as  follows; 

a)  Place  z-transforms  on  a firm  foundation  using  distribu- 
tion theory.  This  has  already  been  done  witli  Laplace  and  Fourier 
transforms  [13].  The  intimate  relation  between  z-transforms  and  the 
Dirac  delta  is  suggested  in  Appendix  A.  Sucli  a foundation  would  allow 
extensions  and  further  development  of  z-transforms, 

b)  The  effects  of  tuning  for  otlier  inputs  and  transfer 
functions  require  analysis.  The  times  when  modified  z-transforms 
and/or  tunable  convolution  are  advantageous  and  in  what  combination 
they  are  advantageous  require  study.  This  would  include  tlie  use 

of  higher  order  holds.  The  first-order  hold  leads  to  Simpson's  rule. 
Equation  (59) . 

Ultimately,  it  is  hoped  that  there  would  be  some  unification 
between  z-transforms,  classical  numerical  methods,  and  distribution 
theory.  And  catastrophe  theory... 


Also,  a literary  production  composed  of  unconnected  parts. 
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Appendix  A.  Z TRANSFORMS  AND  OTHER 
DEFINITIONS  AND  RELATIONSHIPS 


Laplace  transform; 


00 

/ 


-St 


£(f(t))=  / f(t)e  dt 

o 

= f(s) 

Convolution: 


t 

-i 


£(t)  * g(t)s  / f(x)  g(t  - T)dT 
"o 

£ (f(t)  * g(t))  = f(s)  g(s)  . 

The  Dirac  delta,  6(t),  is  the  unity  element  in  convolution. 
Cauchy  Product  of  Power  Series: 

( I 'nV  I S \ = f f . 

\n=0  / \n=0  / n=0  k=0 

Sampling  with  Dirac  delta  distribution: 

00 

f(nT)  = y 6(t  - nT)  f(t)  dt 


z-trans  form; 


Z(£(f(t))) 


00  oo 

-I  / 


6(t  - nT)f(t)e"®*^  dt 


n=0 


-st 


(A-1) 


(A-2) 


(A-3) 


(A-4) 


(A-5) 


(A-6) 


(A-7) 


It  is  noted  that  f(t)e  is  being  sampled,  then  summed.  From  the 
properties  of  the  Dirac  delta. 


45 


/cx> 

b(t  - nT)dc 


(A-8) 


. Z £ (nI)e-“^ 

n=0 

(A-9) 

00 

♦ 

= Z f(nT)z'^ 

n=0 

(A-10) 

= f(z) 

(A-11) 

Equation  (A-10)  seems  to  arise  in  the  literature  like  Venus 

full 

grown  from  the  foam,  or  is  it  like  Athena  from  the  forehead 
Modified  z-transform: 

00 

0®  r 

of  Zeus? 

^ / 5(t  - (n  + q)T)f(t)e"®"  dt  , 0 < q < 1 

n=0  J 
-00 

(A-12) 

: 

00 

' 

f(z)  = Z f(nT  + 
n=0 

(A-13) 

00 

X £<»T  + 
n-K) 

but 

00 

(A-14) 

' 

f(z>  n)  s X + qT)e  ; 

n=0 

(A-15) 

1 

therefore. 

i 

z"  f(z)  = f(z,  q) 

(A-16) 

■1 

where  z*'  is  a fractional  shift. 

* 
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Shifting  theorem: 


00 

! / 

n=0  J 


-St 


6(t  - (n  + m)T)f(t)e  dt  , 0 < m 


Y,  f(nT  + 
n=0 


= Y.  + mT)( 

n=0 


•nsT 


-msT  . 

= e f(z,in) 


= z f(z,m) 
m-1 

= f(z)  - Y,  f(f^T)z"  ; 
n=0 


therefore , 


f(z,m)  = z f(z)  - Y f(t'T)z‘^l 

L n=0  J 


or 


m-1 

z f(z)  = f(z,m)  + z ™ ^ f(nT)z'^ 

n=0 


the  shifting  theorem. 

Shifting  the  other  way. 


00 

f f 

n=0  J 


8(t  - (n  - m)T)f(t)e  dt  , 0 < m 


= X f(nT  - 
n=0 


(A-17) 

(A-18) 

(A-19) 

(A-20) 

(A-21) 

(A-22) 

(A-23) 

(A-24) 


(A-25) 


(A-26) 


1 


i 

j 


j 

I 
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= f(nT  - mT)e“^®^ 

n=0 


= z f(z,-tn) 


-1 

= f(z)  + X 

k=-m 


-ksT 


(A-27) 

(A-28) 

(A-29) 


therefore, 


f ( z , -m)  = z 


f(z)  + X f(kT)z 
k=-m 


(.A- 30) 


or 


z"*  f(z)  = f(z,-m)  - z 2^  z 

k=-m 


-1 

I 


(A-31) 


m-1 

= f(z,-m)  - X = 

n=0 


(A-32) 


Just  because  f(nT),  n < 0,  may  not  be  known  does  not  mean  it  is  zero. 


J 
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Appendix  B.  TRANSFORM  TABLE  FOR  SELECTED  FUNCTIONS* 
AND  DISTRIBUTIONS 


-niealv,  M.,  Tables  of  Laplace.  Heaviside.  Fourier  and  Z-transforms . 
London;  W.  and  R.  Chambers  Ltd.,  1967.  Cadzow,  J.  A.,  Discrete-Time 
Systems . Englewood  Cliffs,  New  Jersey;  Prentice  Hall,  1973. 
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Appendix  C.  DIRAC  DELTA  DISTRIBUTION^ 

f(t)  B(t  - nT)dt  = f(iiT) 


where  f(t)  is  a well-defined  function  at  t = nT. 
S(t)  = 6(-t) 


1 


B(at)  = , , 


&(t) 


5(t  - nT)  ® 

(nl)l  g'(nT)  ^ 0 


t 6(t)  = 0 

f(t)B(t  - nT)  = f(nT)S(t  - nT) 


/ 


B(t  - T)&(t  - nT)dt  = 8(t  - nT) 


00 

/ 


(m) 


(m) 


B(t)  f(t)  = (-1)  f(o) 


00 

/ 


(U 


&(o)  e"®^^  dt  = s 


(C-1) 


(C-2) 

(C-3) 


(C-4) 


(C-5) 

(C-6) 


(C-7) 


(C-8) 


(C-9) 


*Dirac,  P.  A.  M. , The  Principles  of  Quantum  Mechanics,  London; 
Oxford  University  Press,  1958,  Messiah,  A.,  Quantum  Mechanics,  New 
York;  John  Wiley  and  Sons,  1964. 
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Appendix  D.  PAST  VALUES  FOR  DIFFERENCE  EQUATIONS  SOLVING 
DIFFERENTIAL  EQUATIONS  WITH  INITIAL  CONDITIONS 


As  was  seen  in  the  sample  prohUms  in  Siction  , the  order  oi  the 
plant  (and  the  forcing  function)  determiiu-  the  numher  of  start-up  steps 
before  the  recurrence  may  bo  applied.  In  a simulation  containing  differ- 
ent order  differential  equations,  tlie  recurrences  would  start  on  differ- 
ent steps,  which  could  add  complexity  to  the  simulation. 

To  stait  all  recurrences,  regardless  of  their  order  at  Step  1 
(n  = 1),  in  general,  values  before  n = o will  be  required.  First,  the 
required  number  of  steps  must  be  solved  backwards  in  time  before  pio- 
ceeding  forward  in  time  with  the  recurrence  at  Step  1,  The  dilferontial 
equation  could  also  be  solved  backv^;ards  in  time  (odd  derivatives  and 
functions  change  sign,  etc.,  and  the  startup  steps  would  bo  obtained 
for  computing  backwards  as  well  as  the  backwards  recurrence,  which  is 
not  required. 


If  the  forward  start-up  steps  have  already  been  determined,  the 
backwards  ones  may  be  determined  by  substituting  -T  for  T .ind  -n  for  n. 
This  would  also  apply  for  the  recurrence.  For  example,  Kquations  (Ibl), 
(160),  and  (159)  would  become,  respectively. 


X 

-n 


2 X , , . - X + T"  x’  . , - n < -1 

-(n-1)  -(n-2)  -(n-1) 


2 X 1 ■ X - + T^  x’  , , n > 1 
-n+1  -n+2  -n+1 ’ 


(D-1) 


and 


X ••  x , n = o 

o o 


(0-2) 


(D-3) 


With  Equations  (D-2)  and  (D-3),  Equation  (161)  may  then  be  written  as; 


X 

n 


2 X 

n- 


1 


n-2 


+ T“ 


Vi’ 


(D-A) 


Because  the  equations  were  solved  backwards,  the  past  values  required 
for  starting  the  recurrence  at  Step  1 (n  = 1)  are  not  necessarily  the 
past  values  of  the  actual  system.  The  forcing  function  nviy  not  liave 
been  present  before  time  zero  (n  = o) ; nevertheless,  it  must  be  taken 
into  consideration  to  start  the  recurrence  at  Step  1. 
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A nil!..  J example  will  illustrate.  For  x = 1,  x = 1 , x'  = 2. 

o ’ o o ’ 

and  T = I,  Fquations  (159),  (160),  and  (161)  become 

X =1 
o 

Xj^  = l + l'<l+jx2  = 3 

X2=2x3-1  + 1x2  = 7 

Using  Kquations  (D-2),  (f)-3),  and  (D-4)  twice, 
x_j^  = 1-1+|-X2  = l 


X = 1 
o 

x^=2xl-l+lx2=3 
x2=2x3-1  + 1x2  = 7 

The  free  running  oscillator.  Equations  (203),  (204),  and  (205), 
noting  that  cosine  is  an  even  function  and  sine  an  odd  function,  would 
become 

x_^  = x^  (cos  off)  - (sin  uiT)  , n = -1  , (D-5) 

^o  " ’‘o  . n = o , 

and 

x^  = (2  cos  COT)  x^_^  - x^_2,  n > 1 . (U-7) 

In  this  way  a "preprogram”  could  convert  initial  conditions  for  the 
differential  equations  into  the  required  past  values  for  the  difference 
equations  at  Step  1. 

There  is  considerable  confusion  on  this  point  in  the  lituerature 
and  caution  is  advised. 
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